Geoscanning Retailer Database Workflow - VERSION 3

1. Introduction

This document details the code used to intake, clean, geolocate and coalesce administrative data on tobacco retail licenses from three states - Pennsylvania, New Jersey and Delaware. These data will be used to relate the geo-location information provided by GeoScanning project subjects to retailer locations according to experimental parameters to be specified.

This codebase does two things:

1. Creates an initial database of retailers which can be updated as new data are published. (This "baseline" data set consists of 31,034 retailers, of which only 10 could not be geolocated after qualitative research.)

2. Creates a process by which new data sets - published by PA and DE using open source Application Programming Interfaces (APIs), and provided by direct communication from New Jersey - are cross-checked against existing data.

New retailers are added to the database, old ones are marked as closed, repeat entries are logged as continuing their retail activity. All observations are then geo-located and errors checked and corrected - first automatically and then manually. This updating protocol should hypothetically be reproducible.

This document is generally organized by state - with each section documenting a) the creation of a baseline data set and b) the update routine for that state.

The end of the document contains code for making interactive maps and also has a set of data dictionaries.

States can all be updated at once, but more likely they can be done one at a time - code for doing one state in a sitting (uploading the whole baseline set and just replacing one state) is included near the end of this document.

These retailer databases are designed to be integrated into a broader code base which intakes and analyzes geo-location data and relates them to the retailer information. Not only do these data have lat/lon which enable geographic relationships, but they now have a time-window associated with licnesure at that location, allowing for a second check against geo-location info to confirm "exposure" to point-of-sale tobacco locations.

During the summer and fall of 2019, Michael Fichman, city planning researcher and analyst at PennPraxis (of the Weitzman School of Design), took pre-existing codebase reltaed to geo-location data developed for the project and updated it to conform (where possible) to modern tidyverse data wrangling standards and the powerful sf spatial package which is integrated in to the tidyverse and R's powerful visualization package ggplot2.

The next step is to integrate these retailer data, in their now-complete form, to geo-location data, and integrate measures developed by NYU including "stay event" detections of sedentary behavior in addition to previously created mechanisms for spatial joins between retailer buffers and geo-location data.

At the present date, this is a document-in-progress and the code contained heirein is designed to operate on test geo data and a temporary data set of tobacco retailer locations. Subsequent steps to code base development will include the creation of a suite of descriptive statistics for standard data intake and analysis procedures, the creation of an interactive map as part of a subject intake protocol and the development of a "spatial risk surface" analytic approach. Furthermore, retailer information for a multi-state study area will need to be run through the protocol.

1.1. About this code base

This is a "flexible" code base rather than a series of tightly-wrapped functions that take a series of inputs and parameters and output an object. This type of code base is a bit harder to use without a surface understanding of R, but much more adaptable to new intake data types and structures, updates and new analyses. For example, this code will have fewer tightly written functions that can be simply executed, but are very difficult to troubleshoot in the event of bad data or changes to the programming language.

The user needs to be able to operate the R studio environment, execute code in the R studio console, install and load packages, set filepaths for uploading and exporting data, change variable and object names and perform a basic level of troubleshooting using internet forums (like Stack Overflow, github etc) and R documentation if necessary.

That said, this code is open source and subject to deprecation, package updates etc., and is therefore provided with no definitive warranty. This document is designed to annotate the code base to be understandable, reproducible and reflective of research needs.

The elements of this codebase can be used out of order - with states updated one state at a time. Furthermore, the general outline of the leaflet code in can be used at any stage of the process to check the retailer locations for gross errors and explore them more deeply.

1.2. Using This Document

This document is an R Markdown document, which is a publishing language designed embed chunks of R code. It can also be used to create reproducible reports undergirded by R code.

Code can be cut and pasted from this document into a regular R code file or the R console.

Use the table of contents at the left hand side to navigate through this document.

1.3. Workflow

1. Create a baseline data set from previously published license data, standardize the column names and geolocate all of the observations with lat/lon using the Google geocoding API. We refer to this as a "precursor routine" and it is done for each of our three states. This code is shown here for documentation purposes, but shouldn't need to be replicated in the future.

2. Intake new licensing information directly from PA and DE using the Socrata Open Data API (SODA) and from New Jersey through direct communication. The basic workflow here is to read in json data from the API, standardize column names to conform to our standards, convert it to sf data type, timestamp it, write it out with a date stamp for today's date along with the expiration date.

3. Compare the licenses in the new data to those in the baseline data and log their stay in the database. If a licensee is new, create an entry for it. If a licensee no longer has a registration in the new data, mark their license as "EXPIRED". If a licensee is still found and unchanged, extend the time of its stay in the database

4. Geolocate the retailer data by sending any address which lack lat/lon information to the Google geocoding API, which is the industry standard service.

5. Check errors using an automated process - figuring out if there are points which are outside of the state in which they are supposed to be located. Often when there are errors with geolocated data - the errors are quite big - across the country, not just accross the street.

6. Fix errors (bad geocodes and missing address/geo-data) manually using google maps and research using state business databases, business listings and other methods. Periodically create leaflet maps to eye-test the checking.

7. Write out the data for each state to file.

8. Coalesce the data into a multi-state data set and write it to file.

9. Map the data using the leaflet web mapping package for display and use during study intake OR during the error checking process (to examine individual observations.). This leaflet map can be run at any time in order to check errors (see step 6)

2. Setup

2.1. Load Libraries

library(tidyverse)
library(sf)
library(tigris)
library(tidycensus)
library(lubridate)
library(ggmap)
library(jsonlite)
library(RSocrata)
library(readxl)
library(leaflet)
library(leaflet.providers)
library(readxl)

tidyverse - R's mutually intelligible world of data-wrangling and graphics packages including dplyr

sf - The premier spatial package for R, which plays well with ggplot2, dplyr and the other tidyverse packages. Includes geoprocessing capabilities for point and polygon geometries.

tigris - A package for calling shapefiles from the US Dept. of Commerce, including political boundaries, school districts and hydrology

tidycensus - A tidy and sf compatible package for calling the US Census API for any available data from the ACS and decennial census, 2000 to present. Includes geometries from tracts to counties to states if desired.

lubridate - A sophisticated yet simple date-handling package

ggmap - An older mapping graphics package designed for use with ggplot, which has been surpassed by sf but allows for geocoding to APIs such as Google Maps and the Data Science Toolkit.

jsonlite - A package designed to upload, convert and handle JSON files in R

RSocrata - A package used to interact with APIs commonly used to host open data sets. Pennsylvania and Delaware use "SODA" - the Socrata Open Data API as an endpoint for its data.

readxl - A package for reading Excel files to tibble/dataframe

leaflet - One of the most common javascropt/html/css mapping packages for web use. R Markdown supports leaflet functionality and does not require the use of javascript. Leaflet does not suport sf and any vector data beyond points needs to use rgdal, an older spatial package.

leaflet.providers - additional basemaps and styling options for leaflet, including ESRI data

2.2. Register Google API Key

Until recently there were several geocoders of reasonable quality that were available for free, especially for small runs. To deter developers from making endless API calls, high quality services such as the Google geodocidng API are available only for a fee.

This is a geocoding-intensive project, with thousands of observations requiring geocoding. The ggmap package allows for simple geocoding by address in R, with the google API as one of the available extensions.

You can get your own key by signing up for a google developer account, which gets you a 300 dollar starting credit. A thousand requests cost 5 dollars. This should be sufficient.

If you are running this codebase, consult with me for my API key, or get your own.

Information about using the ggmap package and the google API can be found at the 'r-bloggers' site. You can follow that link to get your own account started.

Begin by registering your key with the ggmap function register_google

register_google(key = "YOUR API KEY GOES HERE")

2.3. Create Canonical Names

We create a list of variable names which our final data sets will use. Each state uses different names to describe similar fields (e.g. latitude vs lat) and these standards can change from year to year. We will create our own data standard to which we will force each data set to conform.

The field names are somewhat self-explanatory.

canonical_names <- c("county", "trade_name", "account", "license_type",
                     "expiration_date", "lat", "lon", "address_full",
                     "publish_date", "state", "expired_y_n")

2.4. Load State Shapefiles

In order to double check our geocoded retailer locations to detect serious errors, we load shapefiles for US state boundaries from the US Census Bureau using the tigris package to interact with their API (this does not require a Census API Key). We download the complete shapefile, transform it to an sf object, and then subset each state and reproject the spatial data to WGS84 (aka Web Mercator) using the st_transform function from the sf package.

states_shp <- states()  %>%
  st_as_sf() %>%
  st_transform(crs = 4326)

pa_shp <- states_shp %>%
  filter(STUSPS == "PA")

nj_shp <- states_shp %>%
  filter(STUSPS == "NJ")

de_shp <- states_shp %>%
  filter(STUSPS == "DE")

3. Precursor Routine - Baseline PA Retailer Data

Our baseline PA data come from a shapefile downloaded by the Geoscanning team in 2018 from the Commonwealth's Open Data portal. This is a precursor to our code routine (e.g. our code takes in NEW data). The following steps do not need to be re-run, and act only as documentation for the building of the portion of our retailer database 11/2018-1/2020. Ancestral data are not published by PA, this provides a fairly good look into the recent past.

3.1. Read new data from Socrata

Steps:

1. Grab the data from the PA SODA API - at the time (1/29/2020) the most recent set. (read.socrata)

2. Filter out of state and vending observations

3. Add a time stamp with the date of download (1/29/2020), a column for 'state' and coerce the license expiration date to date ('ymd')

4. Rename relevant columns to lat, lon and address_full

5. Use mutate_if to coerce all factors to character

6. Replace "&" characters with "AND" (better for geocoding)

Address data are in a json-like format (in brackets etc) but this should not bother the Google geocoder.

retailers_Socrata <- read.socrata("https://data.pa.gov/resource/ut72-sft8.json") %>%
  filter(county != "UNKNOWN/OUT OF STATE",
         license_type != "Vending") %>%
  mutate(publish_date = today(tzone = "America/New_York"),
         state = "PA",
         expiration_date = ymd(expiration_date)) %>%
  rename(lat = location_1.latitude,
         lon = location_1.longitude,
         address_full = location_1.human_address) %>%
  dplyr::select(-legal_name, -postal_code, -country) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(address_full = str_replace(address_full, "&", "AND"))

3.2. Read in original shapefile from 2018

The older, shapefile data, are not uniform with the data from the SODA API. We unite the address information, rename various columns to comport with newer data, extract lat/lon values from the geometry, coerce to data frame, remove geometery, convert factors to character (mutate_if).

This code probably won't be necessary to use in the future when you are using downloads from the Socrata SODA API.

retailers_1 <- st_read("//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/inputData/RetailerMaps/RetailerMaps/PA/2018_11/rows_11/rows_11.shp") %>%
  filter(location_3 == "PA",
         county != "UNKNOWN/OUT OF STATE",
         license_ty != "Vending") %>%
  mutate(location_1 = as.character(location_1), 
         location_2 = as.character(location_2), 
         location_3 = as.character(location_3), 
         postal_cod = as.character(postal_cod),
         state = "PA") %>% 
  unite(., "address_full", sep = " ",
        c("location_1", "location_2", "location_3", "postal_cod")) %>%
  dplyr::select(-computed_r, -version, -computed_1, -computed_2, -computed_3, -computed_4,
         -cartodb_id, -country, -location_4, -created_at, -id, -country, -legal_name) %>%
  rename(publish_date = updated_at,
         expiration_date = expiration,
         license_type = license_ty) %>%
  mutate(lon=map_dbl(geometry, ~st_centroid(.x)[[1]]),
         lat=map_dbl(geometry, ~st_centroid(.x)[[2]]))%>%
                  as.data.frame() %>%
                  select(-geometry) %>%
  mutate_if(is.factor, as.character)  %>% 
  mutate(address_full = str_replace(address_full, "&", "AND"))

Check to see how many NA in spatial information by summarizing.

summary(is.na(retailers_Socrata$lat))
summary(is.na(retailers_1$lat))

3.3 Coalesce Data Sets

Check to see which licenses are new in the incoming data. When we find new entries that match entries in our existing data, we update their date of expiry. If data entries in the old database are not represented in the new one, we mark them as expired.

Note the use of the dplyr version of ifelse called if_else here in order to make sure the result is a date and not numeric.

Here's what this does:

1. Full Join of retailers from Nov 2018 (e.g. existing retailers) to retailers_Socrata (e.g. new retailers) joined by both account and trade name. Since the data sets have identical names, this produces an "x" and "y" column for each of our canonical names - x being the left side and y being the right side of the join.

2. If there is no expiration date on file for the Socrata set (e.g. there is no entry for a retailer on file in the newest data set), assume the account is expired. Create a column expired_y_n and mark it "EXPIRED" or else "ACTIVE"

3. If it's expired (e.g. there is no corresponding y entry for a given x entry), mark an expiration date as being the last active date of the license on file from the existing database.

4. Create a new set of canonical columns by checking to see if both x and y columns are filled for each of our variables - if there's no x, put the y data in there and vice versa.

5. Generate lat/lon by checking if either new or existing has missing lat/lon info, if so, take the existing one, if not, take the newest lat/lon.

6. Use similar logic to fill in address, county, publish date and license type data.

7. Remove the columns which are left over from the join (e.g. lat.x, license_type.y) using dplyr::select

In this join the most recent data set goes second (e.g. the "y") data set.

pa_join <- full_join(retailers_1, 
                     retailers_Socrata, by = c("account", "trade_name")) %>%
  mutate(expired_y_n = ifelse(is.na(expiration_date.y) == TRUE, "EXPIRED", "ACTIVE"),
         expiration_date = if_else(expired_y_n == 'EXPIRED', expiration_date.x, ymd(expiration_date.y)),
         lat = ifelse(is.na(lat.x) == FALSE, lat.x, lat.y),
         lon = ifelse(is.na(lon.x) == FALSE, lon.x, lon.y),
         address_full = ifelse(is.na(address_full.x) == FALSE, address_full.x, address_full.y),
         county = ifelse(is.na(county.y) == TRUE, county.x, county.y),
         publish_date = if_else(is.na(publish_date.y) == TRUE, publish_date.x, publish_date.y),
         license_type = ifelse(is.na(license_type.y) == TRUE, license_type.x, license_type.y),
         state = "PA") %>%
  dplyr::select(-lat.x, -lon.x, -lat.y, -lon.y, 
         -address_full.x, -address_full.y, -county.x, -county.y,
         -publish_date.x, -publish_date.y, -state.x, -state.y, -license_type.x, -license_type.y,
         -expiration_date.x, -expiration_date.y)

We again summarize to check on the amount of NA data.

summary(is.na(pa_join$lat))

3.4. Geocode missing information

The google geocoder takes a list of addresses and returns a list of lat/lon pairs.

1. We create a dataframe of observations to be geocoded (to_geocode) which lack lat/lon data.

2. We send that data to the geocoder via ggmap (geocode).

3. We cbind the returned lat/lon pairs to the rest of the columns in the to_geocode dataframe

4. We rbind the rows from our to_geocode set, which now have x-y info, to our data which did not need to be geocoded in the first place and create a dataframe pa_geocoded.

5. Perform a data check to see what might still be missing

to_geocode <- pa_join %>%
  filter(is.na(lat) == TRUE)

pa_geocoded <- geocode(to_geocode$address_full, source =  "google") %>%
  cbind(., to_geocode %>%
          select(-lat, -lon)) %>%
  rbind(pa_join %>%
          filter(is.na(lat) == FALSE), .)

summary(is.na(pa_geocoded$lat))

3.5. Perform quality check for bad geocodes

When geocodes are wrong, they are often very wrong, and set somewhere wildly different from our target area. We do a spatial join with the PA shapefile, keep only the observations which fall outside of PA. We can check these observations and potentially geocode them manually.

We do this through the following steps:

1. Perform a spatial join using st_join from the sf package to see which points from pa_geocoded fall outside PA (excluding those we know are NA at present). To do this we need to reproject our lat/lon data as an sf vector object projected in WGS84 (web mercator, crs = 4326) and join it to pa_shp.

2. We filter those observations which come back as is.na(STUSPS) == TRUE (e.g. they don't join at all to the PA tigris shapefile and don't take on any of it's attributes like STUSPS) and know that these are the ones which fell outside our shapefile.

3. We dplyr::select only our canonical variables from these and extract the lat/lon from their geometries using a mutate statement. We want to turn this back into a dataframe with lat/lon data and not geometries like an sf object. So we set it as.data.frame and remove the geometry column.

4. We now rbind these data back to our pa_geocoded observations which lacked lat/lon in the first place, to create a data set consisting of all our observations which we presume have bad geocodes.

outside_PA <- st_join(pa_geocoded %>% 
          filter(is.na(lat) == FALSE) %>% 
          st_as_sf(., coords = c("lon", "lat"), crs = 4326), 
        pa_shp, 
        join = st_within, 
        left = TRUE) %>%
  filter(is.na(STUSPS) == TRUE) %>%
  dplyr::select(trade_name, account, address_full, county) %>%
  mutate(lon=map_dbl(geometry, ~st_centroid(.x)[[1]]),
         lat=map_dbl(geometry, ~st_centroid(.x)[[2]]))%>%
  as.data.frame() %>%
  dplyr::select(-geometry) %>%
  rbind(., pa_geocoded %>%
          filter(is.na(lat) == TRUE))

3.6. Manually geocode remaining observations

The observations that don't geocode will need to be dealt with manually (e.g. input lat/lon coordinates based on google map research). Similarly, observations that geocode poorly and end up on our outside_PA list need to be manually geocoded.

This has to be done by hand unfortuately. You can go through the observations which either have NA values for lat/lon in pa_geocoded or are found to be outside PA in the outside_PA dataframe. Here we do this with web research and mutate statements where we replace the lat/lon values for individual licenses.

In other data sets and later in this code we export the data to a csv and an RA does the coding, which we then re-import and integrate. The export approach is probably less prone to error from transcription problems.

pa_geocoded_fixed <- pa_geocoded %>%
  mutate(lat = ifelse(account == "51**6778", 39.948989, lat),
         lon = ifelse(account == "51**6778", -75.147568, lon)) %>%
  mutate(lat = ifelse(account == "30**1789", 39.931567, lat),
         lon = ifelse(account == "30**1789", -80.058022, lon)) %>%
  mutate(lat = ifelse(account == "70**2004", 41.341935, lat),
         lon = ifelse(account == "70**2004", -75.992409, lon)) %>%
  mutate(lat = ifelse(account == "51**8445", 39.954983, lat),
         lon = ifelse(account == "51**8445", -75.164616, lon)) %>%
  mutate(lat = ifelse(account == "70**1005", 41.017803, lat),
         lon = ifelse(account == "70**1005", -75.589636, lon)) %>%
  mutate(lat = ifelse(account == "51**3448", 39.950923, lat),
         lon = ifelse(account == "51**3448", -75.163884, lon)) %>%
  mutate(lat = ifelse(account == "64**9445", 41.463524, lat),
         lon = ifelse(account == "64**9445", -75.180809, lon))%>%
  mutate(lat = ifelse(account == "06**1675", 40.557271, lat),
         lon = ifelse(account == "06**1675", -75.977632, lon)) %>%
  mutate(lat = ifelse(account == "09**7593", 40.109250, lat),
         lon = ifelse(account == "09**7593", -74.855736, lon)) %>%
  mutate(lat = ifelse(account == "45**4142", 41.094910, lat),
         lon = ifelse(account == "45**4142", -75.327707, lon)) %>%
  mutate(lat = ifelse(account == "10**7974", 41.071094, lat),
         lon = ifelse(account == "10**7974", -80.045434, lon)) %>%
  mutate(lat = ifelse(account == "17**6746", 40.999937, lat),
         lon = ifelse(account == "17**6746", -78.347428, lon)) %>%
  mutate(lat = ifelse(account == "34**8300", 40.999937, lat),
         lon = ifelse(account == "34**8300", -78.347428, lon)) %>%
  mutate(lat = ifelse(account == "49**1271", 40.687000,  lat), 
         lon = ifelse(account == "49**1271", -76.829610, lon)) %>%
  mutate(lat = ifelse(account == "04**5721", 40.586025,   lat), 
         lon = ifelse(account == "04**5721", -80.257807, lon)) %>%
  mutate(lat = ifelse(account == "51**0499", 39.948303,    lat), 
         lon = ifelse(account == "51**0499", 42.115990,  lon)) %>%
  mutate(lat = ifelse(account == "70**5157", -80.033087,    lat), 
         lon = ifelse(account == "70**5157", -75.164633, lon)) %>% 
  mutate(lat = ifelse(account == "70**1014", 40.960844,    lat), 
         lon = ifelse(account == "70**1014", -77.762074, lon)) %>% 
  mutate(lat = ifelse(account == "70**4002", 41.027503,    lat), 
         lon = ifelse(account == "70**4002", -76.079946, lon)) %>% 
  mutate(lat = ifelse(account == "32**1245", 40.613811, lat), 
         lon = ifelse(account == "32**1245",-79.188575, lon)) %>% 
  mutate(lat = ifelse(account == "51**3963", 39.953296, lat), 
         lon = ifelse(account == "51**3963", -75.171498, lon)) %>%
  mutate(lat = ifelse(account == "45**7969", 41.073454, lat), 
         lon = ifelse(account == "45**7969", -75.572500, lon)) %>% 
  mutate(lat = ifelse(account == "51**6676", 39.947059,  lat), 
         lon = ifelse(account == "51**6676", -75.164530, lon)) %>% 
  mutate(lat = ifelse(account == "14**8220", 40.834775,   lat), 
         lon = ifelse(account == "14**8220", -77.672479, lon)) %>% 
  mutate(lat = ifelse(account == "21**0740", 40.251462, lat), 
         lon = ifelse(account == "21**0740", -76.913070, lon)) %>% 
  mutate(lat = ifelse(account == "70**2075", 39.990787,  lat), 
         lon = ifelse(account == "70**2075", -75.615866, lon)) %>% 
  mutate(lat = ifelse(account == "19**1669", 41.042079,   lat), 
         lon = ifelse(account == "19**1669", -76.304124, lon)) %>% 
  mutate(lat = ifelse(account == "03**4033", 40.813582,    lat), 
         lon = ifelse(account == "03**4033", -79.550406, lon)) %>% 
  mutate(lat = ifelse(account == "70**8042", 41.656431,  lat),
         lon = ifelse(account == "70**8042", -80.424375, lon)) %>%
    mutate(lat = ifelse(account %in% c("70**8010", "70**8009", "70**8006",
                                       "70**8002", "70**8001"), 39.876255,  lat),
         lon = ifelse(account %in% c("70**8010", "70**8009", "70**8006",
                                       "70**8002", "70**8001"), -75.242962, lon)) %>%
    mutate(lat = ifelse(account == "70**7009", 40.282027,  lat),
         lon = ifelse(account == "70**7009", -79.550674, lon)) %>%
  mutate(lat = ifelse(account == "70**5130", 41.995120,   lat),
         lon = ifelse(account == "70**5130", -77.124804, lon)) %>%
  mutate(lat = ifelse(account %in% c("70**4109", "70**4074"), 40.492822, lat),
         lon = ifelse(account %in% c("70**4109", "70**4074"), -77.114086, lon)) %>%
  mutate(lat = ifelse(account == "70**3200", 40.367604,   lat),
         lon = ifelse(account == "70**3200", -78.430881, lon)) %>%
  mutate(lat = ifelse(trade_name == "S & B ASSOCIATES INC", 39.950269,   lat), 
         lon = ifelse(trade_name == "S & B ASSOCIATES INC",  -75.147441, lon)) %>%
  mutate(lat = ifelse(account == "63**9067", 40.159599, lat),
         lon = ifelse(account == "63**9067", -80.276352, lon)) %>%
  mutate(lat = ifelse(account == "63**9165", 40.367822,  lat),
         lon = ifelse(account == "63**9165", -80.290398, lon)) %>%
  mutate(lat = ifelse(account == "23**0772", 39.9049401,  lat),
         lon = ifelse(account == "23**0772", -75.3885344, lon)) %>%
  mutate(lat = ifelse(account == "40**0587", 41.347524,   lat),
         lon = ifelse(account == "40**0587", -75.960558, lon)) %>%
  mutate(lat = ifelse(account == "62**8390", 41.841740,    lat),
         lon = ifelse(account == "62**8390", -79.354088, lon)) %>%
  mutate(lat = ifelse(account == "02**5373", 40.504864,     lat),
         lon = ifelse(account == "02**5373", -79.675642, lon)) %>%
  mutate(lat = ifelse(account == "54**9250", 40.690829, lat),
         lon = ifelse(account == "54**9250", -76.260003, lon)) %>%
  mutate(lat = ifelse(account == "15**9949", 39.986204, lat),
         lon = ifelse(account == "15**9949", -75.805867, lon)) 

# ADD ADDITIONAL GEOCODING MUTATES HERE

3.7. Write out the dataset

We export our PA baseline data set as a csv.

write.csv(pa_geocoded_fixed, "//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/dataOutputs/pa_geocoded_fixed_5_15_20.csv")

4. PA Data - Socrata Routine for New Data

This is the data routine we will use for updating the PA data. Its naming conventions conform to similar routines for NJ and DE and its data outputs will integrate with those other states.

This routine involves downloading license data from the PA SODA API endpoint, (so long as that is the standard format for such data) and cross-checking it against baseline data to close out old retailer locations and add entries for new ones. Observations without a geolocation in the state database are then geocoded and errors are checked.

4.1. Load new data from Socrata

This routine is basically the same as the baseline routine from Section 3.

Steps:

1. Grab the data from the PA SODA API. (read.socrata)

2. Filter to remove out of state and vending observations

3. Add a time stamp with the date of download, a column for 'state' and coerce the license expiration date to date ('ymd')

4. Rename relevant columns to lat, lon and address_full

5. Use mutate_if to coerce all factors to character

6. Replace "&" characters with "AND" (better for geocoding)

Address data are in a json-like format (in brackets etc) but this should not bother the Google geocoder.

retailers_Socrata_pa <- read.socrata("https://data.pa.gov/resource/ut72-sft8.json") %>%
  filter(county != "UNKNOWN/OUT OF STATE",
         license_type != "Vending") %>%
  mutate(publish_date = today(tzone = "America/New_York"),
         state = "PA",
         expiration_date = ymd(expiration_date)) %>%
  rename(lat = location_1.latitude,
         lon = location_1.longitude,
         address_full = location_1.human_address) %>%
  dplyr::select(-legal_name, -postal_code, -country) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(address_full = str_replace(address_full, "&", "AND"))

4.2. Load most current stored data set

We load up our saved baseline data, from wherever that may be stored. When running this code, the filepath in the following code chunk will likely need to be changed to reflect the most recent PA baseline data set.

Note that we dplyr::select only our canonical_names to account for the fact that when we wrote these data out, they likely were imparted with some kind of column denoting a row id. We don't want that.

We also make sure our data types are correct - using a mutate_if statement to turn factors to characters and a mutate statement to turn our dates into date format

stored_pa <- read.csv("//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/dataOutputs/pa_geocoded_fixed_5_15_20.csv") %>%
  dplyr::select(canonical_names) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(expiration_date = ymd(expiration_date),
         publish_date = ymd(publish_date))

4.3. Join new data and stored data

Check to see which licenses are new in the incoming data. When we find new entries that match entries in our existing data, we update their date of expiry. If data entries in the old database are not represented in the new one, we mark them as expired.

Note the use of the dplyr version of ifelse called if_else here in order to make sure the result is a date and not numeric.

Here's what this does:

1. Full Join of retailers from baseline PA data (stored_pa) to new data (retailers_Socrata_pa), joined by both account and trade name. Since the data sets have identical names, this produces an "x" and "y" column for each of our canonical names - x being the left side and y being the right side of the join.

2. If there is no expiration date on file for the new set (e.g. there is no entry for a retailer on file in the newest data set), assume the account is expired - mark it "EXPIRED" or else "ACTIVE"

3. If it's expired (e.g. there is no corresponding y entry for a given x entry), mark an expiration date as being the last active date of the license on file from the existing database.

4. Create a new set of canonical columns by checking to see if both x and y columns are filled for each of our variables - if there's no x, put the y data in there and vice versa.

5. Use similar logic to fill in address, county, publish date and license type data.

6. Remove the columns which are left over from the join (e.g. lat.x, license_type.y) using dplyr::select

In this join the most recent data set goes second (e.g. the "y") data set.

joined_pa <- full_join(stored_pa, 
                     retailers_Socrata_pa, by = c("account", "trade_name")) %>%
  mutate(expired_y_n = ifelse(is.na(expiration_date.y) == TRUE, "EXPIRED", "ACTIVE"),
         expiration_date = if_else(expired_y_n == 'EXPIRED', expiration_date.x, ymd(expiration_date.y)),
         lat = ifelse(is.na(lat.x) == FALSE, lat.x, lat.y),
         lon = ifelse(is.na(lon.x) == FALSE, lon.x, lon.y),
         address_full = ifelse(is.na(address_full.x) == FALSE, address_full.x, address_full.y),
         county = ifelse(is.na(county.y) == TRUE, county.x, county.y),
         publish_date = if_else(is.na(publish_date.y) == TRUE, publish_date.x, publish_date.y),
         license_type = ifelse(is.na(license_type.y) == TRUE, license_type.x, license_type.y),
         state = "PA") %>%
  dplyr::select(-lat.x, -lon.x, -lat.y, -lon.y, 
         -address_full.x, -address_full.y, -county.x, -county.y,
         -publish_date.x, -publish_date.y, -state.x, -state.y, -license_type.x, -license_type.y,
         -expiration_date.x, -expiration_date.y)

Using the summary function, we can see how much NA data we have for our geolocations.

summary(is.na(joined_pa$lat))

4.4. Geocode missing information

This routine is the same as section 3.4, but is re-capitulated here with the object names relevant to this section.

The google geocoder takes a list of addresses and returns a list of lat/lon pairs.

1. We create a dataframe of observations to be geocoded (to_geocode_pa) which lack lat/lon data.

2. We send that data to the geocoder via ggmap (geocode).

3. We cbind the returned lat/lon pairs to the rest of the columns in the to_geocode_pa dataframe

4. We rbind the rows from our to_geocode_pa set, which now have x-y info, to our data which did not need to be geocoded in the first place and create a dataframe geocoded_pa.

5. Perform a data check to see what might still be missing

to_geocode_pa <- joined_pa %>%
  filter(is.na(lat) == TRUE)

geocoded_pa <- geocode(to_geocode_pa$address_full, source =  "google") %>%
  cbind(., to_geocode_pa %>%
          select(-lat, -lon)) %>%
  rbind(joined_pa %>%
          filter(is.na(lat) == FALSE), .)

summary(is.na(geocoded_pa$lat))

4.5. Perform quality check for bad geocodes

When geocodes are wrong, they are often very wrong, and set somewhere wildly different from our target area. We do a spatial join with the PA shapefile, keep only the observations which fall outside of PA. We can check these observations and potentially geocode them manually.

We do this through the following steps:

1. Perform a spatial join using st_join from the sf package to see which points from pa_geocoded fall outside PA (excluding those we know are NA at present). To do this we need to reproject our lat/lon data as an sf vector object projected in WGS84 (web mercator, crs = 4326) and join it to pa_shp.

2. We filter those observations which come back as is.na(STUSPS) == TRUE (e.g. they don't join at all to the PA tigris shapefile and don't take on any of it's attributes like STUSPS) and know that these are the ones which fell outside our shapefile.

3. We dplyr::select only our canonical variables from these and extract the lat/lon from their geometries using a mutate statement. We want to turn this back into a dataframe with lat/lon data and not geometries like an sf object. So we set it as.data.frame and remove the geometry column.

4. We now rbind these data back to our geocoded_pa observations which lacked lat/lon in the first place (filter(is.na(lat) == FALSE)), to recreate our original geocoded_pa data set with a better knowledge of what lacks good lat/lon data. Those missing or bad geocodes are now all labeled NA.

errors_pa <- st_join(geocoded_pa %>% 
          filter(is.na(lat) == FALSE) %>% 
          st_as_sf(., coords = c("lon", "lat"), crs = 4326), 
        pa_shp, 
        join = st_within, 
        left = TRUE) %>%
  filter(is.na(STUSPS) == TRUE) %>%
  mutate(lon=map_dbl(geometry, ~st_centroid(.x)[[1]]),
         lat=map_dbl(geometry, ~st_centroid(.x)[[2]]))%>%
  dplyr::select(canonical_names) %>%
  as.data.frame() %>%
  dplyr::select(-geometry) %>%
  rbind(., geocoded_pa %>%
          filter(is.na(lat) == TRUE))

4.6. Manually geocode remaining observations

The observations that don't geocode will need to be dealt with manually (e.g. input lat/lon coordinates based on google map research, license reasearch etc.,). Similarly, observations that geocode poorly and end up on our errors_pa list need to be manually geocoded.

This has to be done by hand unfortuately.

We write out errors_pa and then type in the lat/lon manually in Excel. This is much easier than doing a series of mutate statements.

The name of this filepath should change each time you do it so as not to overwrite previous work. As you are geocoding the observations by hand, keep the file in csv form, rather than xlsx.

write.csv(errors_pa, "filepath goes here.csv")

Once we've done our geocoding research, and add in the lat/lon data which was missing, we load the (now correct) data back in using a read.csv command and join it to our properly geocoded data. Remember that you need to change the file path from the one in the code chunk below.

We use an anti_join to geocoded to isloate those observations from geocoded_pa that didn't end up our errors_pa data set, and then bind the rows together with the manually corrected errors. Note that inside the rbind function we do some data cleaning with our uploaded errors so that we keep only the canonical names and our column datatypes are in the correct formats.

Note that inside the rbind function we do some data cleaning with our uploaded errors so that we keep only the canonical names and our column datatypes are in the correct formats. There is an important hitch here - when the data were exported to Excel for manual cleaning, they may have had their dates scrambled to a month/day/year format rather than the year/month/day our data use. If an error shows up saying "All formats failed to parse. No formats found." this means that one must adjust the ymd function (which turns a string into a date format) to reflect the correct order (e.g. mdy for month/day/year).

geocoded_pa_fixed <-  geocoded_pa %>% 
  anti_join(., errors_pa) %>%
  rbind(., read.csv("input filepath.csv") %>%
  dplyr::select(canonical_names) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(expiration_date = mdy(expiration_date),
         publish_date = mdy(publish_date)))

4.7. Write out the dataset

This is our final PA data set - geocoded_pa_fixed. It can be written out to file. We will bind it to our DE and NJ data later in the code base.

write.csv(geocoded_pa_fixed, "filepath goes here.csv")

5. Precursor Routine - Baseline NJ Data

The New Jersey data routine is significantly different from the PA and DE routines which use SODA API calls and relatively standardized data. We have NJ data from March of 2018 and August of 2019. The latter data set is thought to be relatively complete. However, New Jersey has not geocoded its data, so this baseline routine both compares the 2018 and 2019 data and integrates them, but it geocodes our baseline data set.

5.1. Read March, 2018 Data

1. We read the 2018 info from excel (read_excel).

2. We pull a five digit zip code out from a nine digit code used in their data set. (see: Zip).

3. We create an address_full field for geocoding and replace "&" characters with "AND".

4. We add columns or change names to correspond to our "canonical" column names.

NJ_3_2018 <- read_excel("//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/inputData/RetailerMaps/RetailerMaps/NJ (ncooper@asc.upenn.edu)/Carreras Tartak J W130294 Provided W130294 - Cigarette retailers.xlsx", skip = 5) %>%
  mutate(Zip = ifelse(str_length(Zip) == 9, str_sub(Zip, end = -5), Zip))%>%
  unite(., "address_full", sep = ", ",
        c("Street", "City", "State", "Zip")) %>%
  mutate(address_full = str_replace(address_full, "&", "AND"))%>%
  mutate(expired_y_n = "ACTIVE", 
         expiration_date = ymd("2019-03-31"),  # is this correct?
         publish_date = ymd("2018-04-01"),
         account = str_c("NJ_18", row_number(), sep = "_"), 
         license_type = NA, county = NA,
         state = "NJ",
         Trade = ifelse(is.na(Trade) == TRUE, Name, Trade)) %>% 
    rename(trade_name = Trade) %>% 
    dplyr::select(-Loc_Code, -Name)

5.2. Geocode Data

We do a preliminary geocode - recall, none of the addresses have lat/lon information, so this procedure is different than the preceding PA sections. We are geocoding all of the data, and we expect to have lots of errors and ambiguities to sort out.

1. We send a vector of addresses to the geocoder via ggmap (geocode).

2. We cbind the returned lat/lon pairs to the rest of the columns in the NJ_3_2018 dataframe, creating a new dataframe NJ_3_2018_geocoded

3. We then check the data to see how many observations failed (summary)

NJ_3_2018_geocoded <- geocode(NJ_3_2018$address_full, source =  "google") %>%
  cbind(., NJ_3_2018)

summary(is.na(NJ_3_2018_geocoded$lat))

5.3. Load August, 2019 Data

Load the second data set using the same routine as in section 5.1

There are almost twice as many observations as our 2018 data. We have come to learn that this is because 2018 is an incomplete data set.

NJ_8_2019 <- read_excel("//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/inputData/RetailerMaps/RetailerMaps/NJ (ncooper@asc.upenn.edu)/Carreras-Tartak J C148764 Cigarette retailers.xlsx", skip = 5)  %>%
  mutate(Zip = ifelse(str_length(Zip) == 9, str_sub(Zip, end = -5), Zip))%>%
  unite(., "address_full", sep = ", ",
        c("Street", "City", "State", "Zip")) %>%
  mutate(address_full = str_replace(address_full, "&", "AND"))%>%
  mutate(expired_y_n = "ACTIVE", 
         expiration_date = ymd("2020-03-31"),
         publish_date = ymd("2019-04-01"),
         account = str_c("NJ_19", row_number(), sep = "_"), 
         license_type = NA, county = NA,
         state = "NJ",
         Trade = ifelse(is.na(Trade) == TRUE, Name, Trade)) %>% 
    rename(trade_name = Trade) %>% 
    dplyr::select(-Loc_Code, -Name)

5.4. Coalesce Data Sets

Join the two data sets, see what's new. If there is a new and old record, assign ID as the oldest of the two account numbers, change the expiration date.

In this join the most recent data set goes second (e.g. the "y") data set.

NJ_join <- full_join(NJ_3_2018_geocoded, NJ_8_2019, by = "address_full") %>%
mutate(expired_y_n = ifelse(is.na(expiration_date.y) == TRUE, "EXPIRED", "ACTIVE"),
         expiration_date = if_else(expired_y_n == 'EXPIRED', expiration_date.x, ymd(expiration_date.y)),
         county = ifelse(is.na(county.y) == TRUE, county.x, county.y),
         publish_date = if_else(is.na(publish_date.y) == TRUE, publish_date.x, publish_date.y),
         license_type = ifelse(is.na(license_type.y) == TRUE, license_type.x, license_type.y),
         account = ifelse(is.na(account.x) == FALSE, account.x, account.y),
         trade_name = ifelse(is.na(trade_name.y) == FALSE, trade_name.y, trade_name.x),
         expiration_date = if_else(is.na(expiration_date.y) == FALSE, expiration_date.y, expiration_date.x),
         expiration_date = if_else(is.na(expiration_date.y) == FALSE, expiration_date.y, expiration_date.x),
         publish_date = if_else(is.na(publish_date.x) == FALSE, publish_date.x, publish_date.y),
         state = "NJ",
         county = NA) %>%
  select(-state.x, - state.y, -county.y, -county.x, -account.x, -account.y, -trade_name.y, -trade_name.x,
         -expired_y_n.x, -expired_y_n.y, -expiration_date.y, -expiration_date.x, -license_type.x, -license_type.y,
         -publish_date.x, -publish_date.y)

We can check our data to see how many unique addresses we have.

NJ_join %>% 
  dplyr::select(address_full) %>% 
  unique() %>% 
  nrow()

NJ_join %>%
  dplyr::select(address_full, trade_name)%>% 
  unique() %>% 
  nrow()

5.5. Geocode missing observations

Now we undertake a second round of geocoding to add lat/lon information to our data set - our newest observations don't have any lat/lon information.

This routine is the same as section 3.4, but is re-capitulated here with the object names relevant to this section.

The google geocoder takes a list of addresses and returns a list of lat/lon pairs.

1. We create a dataframe of observations to be geocoded (to_geocode_nj) which lack lat/lon data.

2. We send that data to the geocoder via ggmap (geocode).

3. We cbind the returned lat/lon pairs to the rest of the columns in the to_geocode_nj dataframe

4. We rbind the rows from our to_geocode_nj set, which now have x-y info, to our data which did not need to be geocoded in the first place and create a dataframe geocoded_nj.

5. Perform a data check to see what might still be missing (summary)

to_geocode_nj <- NJ_join %>%
  filter(is.na(lat) == TRUE)

geocoded_nj <- geocode(to_geocode_nj$address_full, source =  "google") %>%
  cbind(., to_geocode_nj %>%
          select(-lat, -lon)) %>%
  rbind(NJ_join %>%
          filter(is.na(lat) == FALSE), .)

summary(is.na(geocoded_nj$lat))

5.6. Perform data quality check for bad geocodes

This section is roughly the same as 3.5 - the methodology is reprinted here using the appropriate object names

When geocodes are wrong, they are often very wrong, and set somewhere wildly different from our target area. We do a spatial join with the NJ shapefile, keep only the observations which fall outside of NJ. We can check these observations and potentially geocode them manually.

We do this through the following steps:

1. Perform a spatial join using st_join from the sf package to see which points from geocoded_nj fall outside PA (excluding those we know are NA at present). To do this we need to reproject our lat/lon data as an sf vector object projected in WGS84 (web mercator, crs = 4326) and join it to nj_shp.

2. We filter those observations which come back as is.na(STUSPS) == TRUE (e.g. they don't join at all to the NJ tigris shapefile and don't take on any of it's attributes like STUSPS) and know that these are the ones which fell outside our shapefile.

3. We dplyr::select only our canonical variables from these and extract the lat/lon from their geometries using a mutate statement. We want to turn this back into a dataframe with lat/lon data and not geometries like an sf object. So we set it as.data.frame and remove the geometry column.

4. We now rbind these data back to our geocoded_nj observations which lacked lat/lon in the first place, to create a data set consisting of all our observations which we presume have bad geocodes.

errors_nj <- st_join(geocoded_nj %>% 
          filter(is.na(lat) == FALSE) %>% 
          st_as_sf(., coords = c("lon", "lat"), crs = 4326), 
        nj_shp, 
        join = st_within, 
        left = TRUE) %>%
  filter(is.na(STUSPS) == TRUE) %>%
  mutate(lon=map_dbl(geometry, ~st_centroid(.x)[[1]]),
         lat=map_dbl(geometry, ~st_centroid(.x)[[2]]))%>%
  dplyr::select(canonical_names) %>%
  as.data.frame() %>%
  dplyr::select(-geometry) %>%
  rbind(., geocoded_nj %>%
          filter(is.na(lat) == TRUE))

We perform a data check to see how many NA values we still have in our data.

summary(is.na(errors_nj$lat))

nrow(errors_nj)

5.7. Manually geocode missing observations

This routine is similar to Section 4.6 but is re-written here with appropriate object names

The observations that don't geocode will need to be dealt with manually (e.g. input lat/lon coordinates based on google map research, license reasearch etc.,). Similarly, observations that geocode poorly and end up on our errors_nj list need to be manually geocoded.

This has to be done by hand unfortuately.

We write out errors_nj and then type in the lat/lon manually in Excel. This is much easier than doing a series of mutate statements.

The name of this filepath should change each time you do it so as not to overwrite previous work. As you are geocoding the observations by hand, keep the file in csv form, rather than xlsx.

write.csv(errors_nj, "//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/dataOutpurs/errors_nj_6_9_20.csv")

Once we've done our geocoding research, and add in the lat/lon data which was missing, we load the (now correct) data back in using a read.csv command and join it to our properly geocoded data.

We use an anti_join to geocoded to isloate those observations from geocoded_nj that didn't end up our errors_nj data set, and then bind the rows together with the manually corrected errors.

Note that inside the rbind function we do some data cleaning with our uploaded errors so that we keep only the canonical names and our column datatypes are in the correct formats. There is an important hitch here - when the data were exported to Excel for manual cleaning, they may have had their dates scrambled to a month/day/year format rather than the year/month/day our data use. If an error shows up saying "All formats failed to parse. No formats found." this means that one must adjust the ymd function (which turns a string into a date format) to reflect the correct order (e.g. mdy for month/day/year).

We also add a series of mutate statements designed to catch erroneous publish/expiration dates imputed into an old data routine (these data were the ones being used for manual geocoding). This is essentially a series of replacements for old imputed values using ifelse statements to catch and over-write incorrect dates.

Nora Johnston did this manually based on an export of errors_nj, and there are a handful of corrections which I did later where I use mutate commandes to alter the lat/lon info.

geocoded_nj_fixed <-  geocoded_nj %>% 
  anti_join(., errors_nj, by = "account") %>%
  rbind(., read.csv("//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/inputData/RetailerMaps/RetailerMaps/repairedGeocodes/errors_nj_6_9_20_mf_edits.csv") %>%
  dplyr::select(canonical_names) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(expiration_date = mdy(expiration_date),
         publish_date = mdy(publish_date))) %>%
  mutate(publish_date = as.character(publish_date), 
         expiration_date = as.character(expiration_date), 
         publish_date = ifelse(publish_date == '2018-03-08', '2018-04-01', publish_date), 
         publish_date = ifelse(publish_date == '2019-08-01', '2019-04-01', publish_date), 
         expiration_date = ifelse(expiration_date == '2019-03-07','2019-03-31', expiration_date), 
         expiration_date = ifelse(expiration_date == '2020-07-31', '2020-03-31', expiration_date), 
         publish_date = ymd(publish_date), 
         expiration_date = ymd(expiration_date)) %>%
  mutate(lat = ifelse(account == "NJ_19_1794", 40.379592, lat),
         lon = ifelse(account == "NJ_19_1794", -74.580245, lon)) %>%
  mutate(lat = ifelse(account == "NJ_19_797", 40.864672,   lat),
         lon = ifelse(account == "NJ_19_797", -73.961253, lon)) %>%
  mutate(lat = ifelse(account == "NJ_18_2339", 40.741128,   lat),
         lon = ifelse(account == "NJ_18_2339", -74.535828, lon))

5.8. Write out the data

This is our final NJ data set - geocoded_nj_fixed. It can be written out to file. We will bind it to our DE and PA data later in the code base.

write.csv(geocoded_nj_fixed, "//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/dataOutputs/nj_geocoded_7_20_20.csv")

6. Routine for New NJ Data (Expected August, 2020)

Presumably, New Jersey will be releasing new data in August, 2020. Because their data sets have been somewhat non-standard to this point, what follows is a possible workflow to intake and clean the new NJ data, presuming the standards are exactly the same of last year.

6.1. Load Baseline data

This is our standard loading procedure, which keeps only the relevant column names and makes sure dates and character formats are correct. If this creates an error, the data may have been disturbed.

stored_nj <- read.csv("filepath goes here") %>%
  dplyr::select(canonical_names) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(expiration_date = ymd(expiration_date),
         publish_date = ymd(publish_date))

6.2. Load New Data

We load new data presuming same data specifications as August 2019 data (same columns and data types.).

This is a spot where there is potential trouble.

One will need to specify the filepath and in the event these excel data are not formatted the same as last year, alter the skip command in the read_excel function. skip specifies how many rows to skip over before finding the column headers. Look over the data set in excel to see how it looks.

If column names are different from the assumed format, using rename commands in this function or changing column names in the excel sheet should solve the problem.

retailers_new_NJ <- read_excel("//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/inputData/RetailerMaps/RetailerMaps/NJ (ncooper@asc.upenn.edu)/NJ Cigarette Retailers 2020 (OPRA C162332) (1).xlsx", skip = 5)  %>%
  mutate(Zip = ifelse(str_length(Zip) == 9, str_sub(Zip, end = -5), Zip))%>%
  unite(., "address_full", sep = ", ",
        c("Street", "City", "State", "Zip")) %>%
  mutate(address_full = str_replace(address_full, "&", "AND"))%>%
  mutate(expired_y_n = "ACTIVE", 
         expiration_date = ymd("2021-03-31"),  # is this correct?
         publish_date = ymd("2020-04-01"),
         account = str_c("NJ_20", row_number(), sep = "_"), 
         license_type = NA, county = NA,
         state = "NJ",
         Trade = ifelse(is.na(Trade) == TRUE, Name, Trade)) %>% 
    rename(trade_name = Trade) %>% 
    dplyr::select(-Loc_Code, -Name)

6.3. Coalesce Data Sets

Check to see which licenses are new. If they are old, we can throw them out and when we add the new data, we will just update their date of expiry. If data in the old database are not represented in the new one, we keep them in, but mark them as expired.

Note the use of the dplyr version of ifelse called if_else here in order to make sure the result is a date and not numeric.

Here's what this does:

1. Full Join of retailers from baseline NJ data (stored_nj) to new data (retailers_new_NJ), joined by both account and trade name. Since the data sets have identical names, this produces an "x" and "y" column for each of our canonical names - x being the left side and y being the right side of the join.

2. If there is no expiration date on file for the new set (e.g. there is no entry for a retailer on file in the newest data set), assume the account is expired - mark it "EXPIRED" or else "ACTIVE"

3. If it's expired (e.g. there is no corresponding y entry for a given x entry), mark an expiration date as being the last active date of the license on file from the existing database.

4. Create a new set of canonical columns by checking to see if both x and y columns are filled for each of our variables - if there's no x, put the y data in there and vice versa.

5. Use similar logic to fill in address, county, publish date and license type data.

6. Remove the columns which are left over from the join (e.g. lat.x, license_type.y) using dplyr::select

In this join the most recent data set goes second (e.g. the "y") data set.

NJ_join <- full_join(stored_nj, retailers_new_NJ, by = "address_full") %>%
mutate(expired_y_n = ifelse(is.na(expiration_date.y) == TRUE, "EXPIRED", "ACTIVE"),
         expiration_date = if_else(expired_y_n == 'EXPIRED', expiration_date.x, ymd(expiration_date.y)),
         county = ifelse(is.na(county.y) == TRUE, county.x, county.y),
         publish_date = if_else(is.na(publish_date.y) == TRUE, publish_date.x, publish_date.y),
         license_type = ifelse(is.na(license_type.y) == TRUE, license_type.x, license_type.y),
         account = ifelse(is.na(account.x) == FALSE, account.x, account.y),
         trade_name = ifelse(is.na(trade_name.y) == FALSE, trade_name.y, trade_name.x),
         expiration_date = if_else(is.na(expiration_date.y) == FALSE, expiration_date.y, expiration_date.x),
         expiration_date = if_else(is.na(expiration_date.y) == FALSE, expiration_date.y, expiration_date.x),
         publish_date = if_else(is.na(publish_date.x) == FALSE, publish_date.x, publish_date.y),
         state = "NJ",
         county = NA) %>%
  select(-state.x, - state.y, -county.y, -county.x, -account.x, -account.y, -trade_name.y, -trade_name.x,
         -expired_y_n.x, -expired_y_n.y, -expiration_date.y, -expiration_date.x, -license_type.x, -license_type.y,
         -publish_date.x, -publish_date.y) %>%
  dplyr::select(canonical_names)

6.4. Geocode missing observations

Now we geocode the observations which lack a lat/lon (to_geocode_nj) by sending addresses to the Google geocoding engine, fetching a lat/lon pair (geocode) and then binding the lat/lon columns back to the rest of the data (cbind) and binding the geocoded observations back to those data which had a lat/lon value (rbind).

The summary function will tell us how many geocodes have been unsuccessful up until this point.

to_geocode_nj <- NJ_join %>%
  filter(is.na(lat) == TRUE)

geocoded_nj <- geocode(to_geocode_nj$address_full, source =  "google") %>%
  cbind(., to_geocode_nj %>%
          select(-lat, -lon)) %>%
  rbind(NJ_join %>%
          filter(is.na(lat) == FALSE), .)

summary(is.na(geocoded_nj$lat))

6.5. Perform data quality check for bad geocodes

This section is roughly the same as 3.5 - the methodology is reprinted here using the appropriate object names

When geocodes are wrong, they are often very wrong, and set somewhere wildly different from our target area. We do a spatial join with the NJ shapefile, keep only the observations which fall outside of NJ. We can check these observations and potentially geocode them manually.

We do this through the following steps:

1. Perform a spatial join using st_join from the sf package to see which points from geocoded_nj fall outside NJ (excluding those we know are NA at present). To do this we need to reproject our lat/lon data as an sf vector object projected in WGS84 (web mercator, crs = 4326) and join it to nj_shp.

2. We filter those observations which come back as is.na(STUSPS) == TRUE (e.g. they don't join at all to the NJ tigris shapefile and don't take on any of it's attributes like STUSPS) and know that these are the ones which fell outside our shapefile.

3. We dplyr::select only our canonical variables from these and extract the lat/lon from their geometries using a mutate statement. We want to turn this back into a dataframe with lat/lon data and not geometries like an sf object. So we set it as.data.frame and remove the geometry column.

4. We now rbind these data back to our geocoded_nj observations which lacked lat/lon in the first place, to create a data set consisting of all our observations which we presume have bad geocodes.

errors_nj <- st_join(geocoded_nj %>% 
          filter(is.na(lat) == FALSE) %>% 
          st_as_sf(., coords = c("lon", "lat"), crs = 4326), 
        nj_shp, 
        join = st_within, 
        left = TRUE) %>%
  filter(is.na(STUSPS) == TRUE) %>%
  mutate(lon=map_dbl(geometry, ~st_centroid(.x)[[1]]),
         lat=map_dbl(geometry, ~st_centroid(.x)[[2]]))%>%
  dplyr::select(canonical_names) %>%
  as.data.frame() %>%
  dplyr::select(-geometry) %>%
  rbind(., geocoded_nj %>%
          filter(is.na(lat) == TRUE)%>%
          dplyr::select(canonical_names))

We perform a data check to see how many NA observations remain in our data.

summary(is.na(errors_nj$lat))

nrow(errors_nj)

6.6. Manually geocode missing observations

This routine is similar to Section 4.6 but is re-written here with appropriate object names

The observations that don't geocode will need to be dealt with manually (e.g. input lat/lon coordinates based on google map research, license reasearch etc.,). Similarly, observations that geocode poorly and end up on our errors_nj list need to be manually geocoded.

This has to be done by hand unfortuately.

We write out errors_nj and then type in the lat/lon manually in Excel. This is much easier than doing a series of mutate statements.

The name of this filepath should change each time you do it so as not to overwrite previous work. As you are geocoding the observations by hand, keep the file in csv form, rather than xlsx.

write.csv(errors_nj, "filepath goes here")

Once we've done our geocoding research, and add in the lat/lon data which was missing, we load the (now correct) data back in using a read.csv command and join it to our properly geocoded data.

We use an anti_join to geocoded to isloate those observations from geocoded_nj that didn't end up our errors_nj data set, and then bind the rows together with the manually corrected errors. Note that inside the rbind function we do some data cleaning with our uploaded errors so that we keep only the canonical names and our column datatypes are in the correct formats.

geocoded_nj_fixed <-  geocoded_nj %>% 
  anti_join(., errors_nj) %>%
  rbind(., read.csv("FILEPATH FOR HAND-GEOCODED OBSERVATIONS") %>%
  dplyr::select(canonical_names) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(expiration_date = ymd(expiration_date),
         publish_date = ymd(publish_date)))

6.7. Write out the data

This is our final NJ data set - geocoded_nj_fixed. It can be written out to file. We will bind it to our DE and PA data later in the code base.

write.csv(geocoded_nj_fixed, "filepath goes here")

7. Precursor Routine - Baseline Delaware data

This precursor routine was used to load and clean the first set of Delaware data we collected - in March of 2020. This is a precursor to our code routine (e.g. our code takes in NEW data). The following steps do not need to be re-run, and act only as documentation for the building of the portion of our retailer database.

Comes from here: https://data.delaware.gov/Licenses-and-Certifications/Delaware-Business-Licenses/5zy2-grhr

There are many license types in these data, and those tobacco related include wholesalers and internet retailers

This will need to be filtered (by category) on the instructions of the team.

We should also filter out wholsalers - many of these are outside delaware - I suppose need to be registered there.

7.1. Load Data From Socrata

Data from March, 2020

This routine is basically the same as the baseline routine from PA, which uses the same general data standards, but it applies to a broad set of business licenses that are published as a group.

Steps:

1. Grab the data from the DE SODA API. (read.socrata)

2. Filter to keep only data of the categories TOBACCO PRODUCTS RETAILER and RETAILER-TOBACCO

3. Add a time stamp with the date of download, a column for 'state' and coerce the license expiration date to date ('ymd')

4. Rename relevant columns to lat, lon, account and license_type

5. unite address fields into one address column for geocoding

6. Remove un-needed variables using select

7. Use mutate_if to coerce all factors to character

8. Replace "&" characters with "AND" (better for geocoding)

retailers_Socrata_de_pre <- read.socrata("https://data.delaware.gov/resource/5zy2-grhr.json") %>%
  filter(str_detect(category, "TOBACCO PRODUCTS RETAILER") | str_detect(category, "RETAILER-TOBACCO")) %>%
  mutate(publish_date = ymd(current_license_valid_from),
         expiration_date = ymd(current_license_valid_to),
         expired_y_n = "ACTIVE",
         trade_name = business_name,
         county = NA,
         state = "DE") %>%
  rename(lat = geocoded_location.latitude,
         lon = geocoded_location.longitude,
         account = license_number,
         license_type = category) %>%
  mutate(zip = ifelse(str_length(zip) == 9, str_sub(zip, end = -5), zip)) %>%
  unite(., "address_full", sep = " ",
        c("address_2", "city", "state", "zip")) %>%
  mutate(state = "DE") %>%
  dplyr::select(-business_name, -current_license_valid_from, 
                -current_license_valid_to, -country, -geocoded_location.human_address,
                -address_1) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(address_full = str_replace(address_full, "&", "AND"))

Use the summary command to explore NA data

summary(is.na(retailers_Socrata_de_pre$lat))

7.2. Geocode data

Now we geocode the observations which lack a lat/lon (to_geocode_de_pre) by sending addresses to the Google geocoding engine, fetching a lat/lon pair (geocode) and then binding the lat/lon columns back to the rest of the data (cbind) and binding the geocoded observations back to those data which had a lat/lon value (rbind).

The summary function will tell us how many geocodes have been unsuccessful up until this point.

to_geocode_de_pre <- retailers_Socrata_de_pre %>%
  filter(is.na(lat) == TRUE)

de_geocoded_pre <- geocode(to_geocode_de_pre$address_full, source =  "google") %>%
  cbind(., to_geocode_de %>%
          dplyr::select(-lat, -lon)) %>%
  rbind(retailers_Socrata_de_pre %>%
          filter(is.na(lat) == FALSE), .)

summary(is.na(de_geocoded_pre$lat))

7.3. Perform data quality check for bad geocodes

One strange issue with DE data is that there are observations in there that are absolutely outside Delaware - the licenses themselves say so. We clean these out later, but they obviously trip up the following routine.

When geocodes are wrong, they are often very wrong, and set somewhere wildly different from our target area. We do a spatial join with the DE shapefile, keep only the observations which fall outside of DE. We can check these observations and potentially geocode them manually.

We do this through the following steps:

1. Perform a spatial join using st_join from the sf package to see which points from de_geocoded_pre fall outside DE (excluding those we know are NA at present). To do this we need to reproject our lat/lon data as an sf vector object projected in WGS84 (web mercator, crs = 4326) and join it to de_shp.

2. We filter those observations which come back as is.na(STUSPS) == TRUE (e.g. they don't join at all to the DE tigris shapefile and don't take on any of it's attributes like STUSPS) and know that these are the ones which fell outside our shapefile.

3. We dplyr::select only our canonical variables from these and extract the lat/lon from their geometries using a mutate statement. We want to turn this back into a dataframe with lat/lon data and not geometries like an sf object. So we set it as.data.frame and remove the geometry column.

4. We now rbind these data back to our de_geocoded_pre observations which lacked lat/lon in the first place, to create a data set consisting of all our observations which we presume have bad geocodes.

there are LOTS of addresses which are legitimately outside delaware (check the zip codes, they are in other states)

errors_de_pre <- st_join(de_geocoded_pre %>% 
          filter(is.na(lat) == FALSE) %>% 
          st_as_sf(., coords = c("lon", "lat"), crs = 4326), 
        de_shp, 
        join = st_within, 
        left = TRUE) %>%
  filter(is.na(STUSPS) == TRUE) %>%
  mutate(lon=map_dbl(geometry, ~st_centroid(.x)[[1]]),
         lat=map_dbl(geometry, ~st_centroid(.x)[[2]]))%>%
  dplyr::select(canonical_names) %>%
  as.data.frame() %>%
  dplyr::select(-geometry) %>%
  rbind(., de_geocoded_pre %>%
          filter(is.na(lat) == TRUE))

We check our data to see how many NA values came back.

summary(is.na(errors_de_pre$lat))

7.4. Manually geocode remaining observations

The observations that don't geocode will need to be dealt with manually (e.g. input lat/lon coordinates based on google map research). Similarly, observations that geocode poorly and end up on our errors_de_pre list need to be manually geocoded.

This has to be done by hand unfortuately. You can go through the observations which either have NA values for lat/lon in de_geocoded_pre or are found to be outside DE in the errors_de_pre dataframe. Here we do this with web research and mutate statements where we replace the lat/lon values for individual licenses.

In other data sets and later in this code we export the data to a csv and an RA does the coding, which we then re-import and integrate. The export approach is probably less prone to error from transcription problems.

There are ABSOLUTELY going to be licenses from outside Delaware on this list. You have to check the zip codes.

Note that in our last step, we remove any obsrvations that are outside the state with a filter command.

de_geocoded_pre_fixed <- de_geocoded_pre %>%
  mutate(lat = ifelse(account == 2004200459, 39.563789,  lat),
         lon = ifelse(account == 2004200459, -75.653397, lon))  %>%
  mutate(lat = ifelse(account == 2000107842, 38.925108,   lat),
         lon = ifelse(account == 2000107842, -75.567057, lon)) %>%
  mutate(lat = ifelse(account == 2016604575, 39.649600,   lat),
         lon = ifelse(account == 2016604575, -75.621843 , lon)) %>%
  mutate(lat = ifelse(account == 2007220039, 39.739074,   lat),
         lon = ifelse(account == 2007220039, -75.623293 , lon)) %>%
  mutate(lat = ifelse(account == 1989005158, 39.823859,   lat),
         lon = ifelse(account == 1989005158, -75.499926 , lon)) %>%
  mutate(lat = ifelse(account == 2019603108, 38.455972,   lat),
         lon = ifelse(account == 2019603108, -75.559722 , lon)) %>%
  filter(!account %in% c(2016605565, 2012101218, 2010102260, 1998204762, 2019100273,
                         1989023113, 2020701642, 2020701644, 2018602428, 2018600730,
                         2016100150, 1996105149))

7.5. Write out the dataset

We export our DE baseline data set as a csv.

write.csv(de_geocoded_pre_fixed, "//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/dataOutputs/de_geocoded_pre_6_8_20.csv")

8. DE Data - Socrata Routine for New Data

This is the data routine we will use for updating the DE data. Its naming conventions conform to similar routines for NJ and PA and its data outputs will integrate with those other states.

This routine involves downloading license data from the DE SODA API endpoint, (so long as that is the standard format for such data) and cross-checking it against baseline data to close out old retailer locations and add entries for new ones. Observations without a geolocation in the state database are then geocoded and errors are checked.

8.1. Load new data from Socrata

This routine is basically the same as the baseline routine from PA, which uses the same general data standards, but it applies to a broad set of business licenses that are published as a group.

Steps:

1. Grab the data from the DE SODA API. (read.socrata)

2. Filter to keep only data of the categories TOBACCO PRODUCTS RETAILER and RETAILER-TOBACCO

3. Add a time stamp with the date of download, a column for 'state' and coerce the license expiration date to date ('ymd')

4. Rename relevant columns to lat, lon, account and license_type

5. Filter just for observations with zip codes between 197001 and 19980 (Delaware Zip codes) - this avoids problems we had creating baseline data with odd out-of-state points

6. unite address fields into one address column for geocoding

7. Remove un-needed variables using select

8. Use mutate_if to coerce all factors to character

9. Replace "&" characters with "AND" (better for geocoding)

retailers_Socrata_de <- read.socrata("https://data.delaware.gov/resource/5zy2-grhr.json") %>%
  filter(str_detect(category, "TOBACCO PRODUCTS RETAILER") | str_detect(category, "RETAILER-TOBACCO")) %>%
  mutate(publish_date = ymd(current_license_valid_from),
         expiration_date = ymd(current_license_valid_to),
         expired_y_n = "ACTIVE",
         trade_name = business_name,
         county = NA,
         state = "DE") %>%
  rename(lat = geocoded_location.latitude,
         lon = geocoded_location.longitude,
         account = license_number,
         license_type = category) %>%
  mutate(zip = ifelse(str_length(zip) == 9, str_sub(zip, end = -5), zip)) %>%
  filter(zip >= 19701 & zip <= 19980) %>%
  unite(., "address_full", sep = " ",
        c("address_2", "city", "state", "zip")) %>%
  mutate(state = "DE") %>%
  dplyr::select(-business_name, -current_license_valid_from, 
                -current_license_valid_to, -country, -geocoded_location.human_address,
                -address_1) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(address_full = str_replace(address_full, "&", "AND"))

8.2. Load most current stored data set

We load up our saved baseline data, from wherever that may be stored. When running this code, the filepath in the following code chunk will likely need to be changed to reflect the most recent DE baseline data set.

Note that we dplyr::select only our canonical_names to account for the fact that when we wrote these data out, they likely were imparted with some kind of column denoting a row id. We don't want that.

We also make sure our data types are correct - using a mutate_if statement to turn factors to characters and a mutate statement to turn our dates into date format

stored_de <- read.csv("//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/dataOutputs/de_geocoded_pre_6_8_20.csv") %>%
  dplyr::select(canonical_names) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(account = as.character(account)) %>%
  mutate(expiration_date = ymd(expiration_date),
         publish_date = ymd(publish_date))

8.3. Join new data and stored data

Check to see which licenses are new in the incoming data. When we find new entries that match entries in our existing data, we update their date of expiry. If data entries in the old database are not represented in the new one, we mark them as expired.

Note the use of the dplyr version of ifelse called if_else here in order to make sure the result is a date and not numeric.

Here's what this does:

1. Full Join of retailers from baseline DE data (stored_de) to new data (retailers_Socrata_de), joined by both account and trade name. Since the data sets have identical names, this produces an "x" and "y" column for each of our canonical names - x being the left side and y being the right side of the join.

2. If there is no expiration date on file for the new set (e.g. there is no entry for a retailer on file in the newest data set), assume the account is expired - mark it "EXPIRED" or else "ACTIVE"

3. If it's expired (e.g. there is no corresponding y entry for a given x entry), mark an expiration date as being the last active date of the license on file from the existing database.

4. Create a new set of canonical columns by checking to see if both x and y columns are filled for each of our variables - if there's no x, put the y data in there and vice versa.

5. Use similar logic to fill in address, county, publish date and license type data.

6. Remove the columns which are left over from the join (e.g. lat.x, license_type.y) using dplyr::select

In this join the most recent data set goes second (e.g. the "y") data set.

joined_de <- full_join(stored_de, 
                     retailers_Socrata_de, by = c("account", "trade_name")) %>%
  mutate(expired_y_n = ifelse(is.na(expiration_date.y) == TRUE, "EXPIRED", "ACTIVE"),
         expiration_date = if_else(expired_y_n == 'EXPIRED', expiration_date.x, ymd(expiration_date.y)),
         lat = ifelse(is.na(lat.x) == FALSE, lat.x, lat.y),
         lon = ifelse(is.na(lon.x) == FALSE, lon.x, lon.y),
         address_full = ifelse(is.na(address_full.x) == FALSE, address_full.x, address_full.y),
         county = ifelse(is.na(county.y) == TRUE, county.x, county.y),
         publish_date = if_else(is.na(publish_date.y) == TRUE, publish_date.x, publish_date.y),
         license_type = ifelse(is.na(license_type.y) == TRUE, license_type.x, license_type.y),
         state = "DE") %>%
  dplyr::select(-lat.x, -lon.x, -lat.y, -lon.y, 
         -address_full.x, -address_full.y, -county.x, -county.y, -expired_y_n.x, -expired_y_n.y,
         -publish_date.x, -publish_date.y, -state.x, -state.y, -license_type.x, -license_type.y,
         -expiration_date.x, -expiration_date.y)

A summary call determines if there are any bad geocodes.

summary(is.na(joined_de$lat))

If there are NO TRUE values to this statement, the following step can be skipped (although you will need to rename joined_de as geocoded_de to get step 8.5 to work correctly).

8.4. Geocode missing information

Now we geocode the observations which lack a lat/lon (to_geocode_de) by sending addresses to the Google geocoding engine, fetching a lat/lon pair (geocode) and then binding the lat/lon columns back to the rest of the data (cbind) and binding the geocoded observations back to those data which had a lat/lon value (rbind).

The summary function will tell us how many geocodes have been unsuccessful up until this point.

to_geocode_de <- joined_de %>%
  filter(is.na(lat) == TRUE)

geocoded_de <- geocode(to_geocode_de$address_full, source =  "google") %>%
  cbind(., to_geocode_de %>%
          select(-lat, -lon)) %>%
  rbind(joined_de %>%
          filter(is.na(lat) == FALSE), .)

summary(is.na(geocoded_de$lat))

If this step gives you an error reading Error in if (location == "") return(return_failed_geocode(output)) : argument is of length zero, run the following code:

geocoded_de <- joined_de

8.5. Perform quality check for bad geocodes

When geocodes are wrong, they are often very wrong, and set somewhere wildly different from our target area. We do a spatial join with the DE shapefile, keep only the observations which fall outside of DE. We can check these observations and potentially geocode them manually.

We do this through the following steps:

1. Perform a spatial join using st_join from the sf package to see which points from geocoded_de fall outside DE (excluding those we know are NA at present). To do this we need to reproject our lat/lon data as an sf vector object projected in WGS84 (web mercator, crs = 4326) and join it to de_shp.

2. We filter those observations which come back as is.na(STUSPS) == TRUE (e.g. they don't join at all to the DE tigris shapefile and don't take on any of it's attributes like STUSPS) and know that these are the ones which fell outside our shapefile.

3. We dplyr::select only our canonical variables from these and extract the lat/lon from their geometries using a mutate statement. We want to turn this back into a dataframe with lat/lon data and not geometries like an sf object. So we set it as.data.frame and remove the geometry column.

4. We now rbind these data back to our geocoded_de observations which lacked lat/lon in the first place, to create a data set consisting of all our observations which we presume have bad geocodes.

errors_de <- st_join(geocoded_de %>% 
          filter(is.na(lat) == FALSE) %>% 
          st_as_sf(., coords = c("lon", "lat"), crs = 4326), 
        de_shp, 
        join = st_within, 
        left = TRUE) %>%
  filter(is.na(STUSPS) == TRUE) %>%
  mutate(lon=map_dbl(geometry, ~st_centroid(.x)[[1]]),
         lat=map_dbl(geometry, ~st_centroid(.x)[[2]]))%>%
  dplyr::select(canonical_names) %>%
  as.data.frame() %>%
  dplyr::select(-geometry) %>%
  rbind(., geocoded_de %>%
          filter(is.na(lat) == TRUE))

8.6. Manually geocode remaining observations

The observations that don't geocode will need to be dealt with manually (e.g. input lat/lon coordinates based on google map research). Similarly, observations that geocode poorly and end up on our errors_de_pre list need to be manually geocoded.

This has to be done by hand unfortuately. You can go through the observations which either have NA values for lat/lon in de_geocoded_pre or are found to be outside DE in the errors_de_pre dataframe. Here we do this with web research and mutate statements where we replace the lat/lon values for individual licenses.

In other data sets and elsewhere in this code we export the data to a csv and an RA does the coding, which we then re-import and integrate. The export approach is probably less prone to error from transcription problems.

geocoded_de_fixed <- geocoded_de %>%
  mutate(lat = ifelse(account == "2003101733", 38.816784,   lat),
         lon = ifelse(account == "2003101733", -75.299722, lon)) %>%
  mutate(lat = ifelse(account == "2020703345", 38.456186,     lat),
         lon = ifelse(account == "2020703345", -75.573836, lon)) %>%
  mutate(lat = ifelse(account == "2010104755", 38.456700,     lat),
         lon = ifelse(account == "2010104755", -75.555881, lon)) %>%
  mutate(lat = ifelse(account == "2011119693", 38.455830,    lat),
         lon = ifelse(account == "2011119693", -75.559666, lon)) %>%
  mutate(lat = ifelse(account == "2000101625", 38.458207,    lat),
         lon = ifelse(account == "2000101625", -75.559746, lon)) %>%
  mutate(lat = ifelse(account == "1989033057", 39.070951,     lat),
         lon = ifelse(account == "1989033057", -75.554093, lon))  %>%
  mutate(lat = ifelse(account == "2019101246", 39.218640,    lat),
         lon = ifelse(account == "2019101246", -75.571138, lon))%>%
  mutate(lat = ifelse(account == "2013103041", 38.509547,     lat),
         lon = ifelse(account == "2013103041", -75.237918, lon))       

8.7. Write out the dataset

This is our final DE data set - geocoded_de_fixed. It can be written out to file. We will bind it to our NJ and PA data later in the code base.

If you are running this code base to do an update, change the filepath below

write.csv(geocoded_de_fixed, "//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/dataOutputs/de_geocoded_10_13_20.csv")

9. Coalescing the master data set

If we are running a routine for all of the states at once, we can coalesce the three states and write out the data set.

allStates <- rbind(geocoded_de_fixed, geocoded_nj_fixed) %>%
  rbind(., geocoded_pa_fixed)

write.csv(allStates, "//jove.design.upenn.edu/Dept-Shares/prax/01 Project Folders/2019_Annenberg_GeoScanning/dataOutputs/all_Retailers_10_20_20.csv")

9.1 Methodology for Updating One State at a Time

It's likely that one would just run a routine for a single state and then update a portion of the data set.

The methodology for doing this would be as follows:

1. Run the routine for the state of interest. (The below example uses NJ)

2. Upload the allStates data set (e.g. the most current, multi-state set) and apply the cleaning convention so that the data types are correct.

3. Add the new observations by filtering out the observations from the state of interest from allStates and using an rbind command to append the new NJ data set.

4. Write the data to file.

allStates_updated <- read.csv("filepath goes here") %>%
  dplyr::select(canonical_names) %>%
  mutate_if(is.factor, as.character) %>%
  mutate(account = as.character(account)) %>%
  mutate(expiration_date = ymd(expiration_date),
         publish_date = ymd(publish_date)) %>%
  filter(state != "NJ") %>%
  rbind(., geocoded_nj_fixed)

write.csv(allStates_updated, "filepath goes here")

10. Make A Leaflet Map

The data can be visualized using the leaflet webmapping package for R Markdown. The following map subdivides the data into active and expired licenses (selectable using a menu at the top right) and allows the user to hover over any point and view the name and address of the license holder. The basemaps are open source and available in the leaflet.providers library.

Additional features can be added to the map, including but not limited to search functions, polygon geometries, and more.

R-Markdown documents can be knit easily into portable html files which can be sent and received easily due to their relatively small file size.

For elaborate data dashboarding, R-Shiny provides more robust joining and selecting functionality, integrated with active charts and tables.

active <- allStates %>%
  filter(expired_y_n == "ACTIVE")

expired <- allStates %>%
  filter(expired_y_n != "ACTIVE")

m4 <- leaflet() %>%
  setView(lng = -77.475127, lat = 40.726604, zoom = 07) %>%
  addProviderTiles(providers$Esri.WorldTopoMap) %>%
  addCircleMarkers(lng= as.numeric(active$lon), 
                   lat= as.numeric(active$lat), 
                   radius = 1,
                   fillOpacity = 1,
                   color= "red",
                   label = paste(active$trade_name, "|",
                           active$address_full),
                           group = "Active Retailers") %>%
  addCircleMarkers(lng= as.numeric(expired$lon), 
                   lat= as.numeric(expired$lat), 
                   radius = 0.5,
                   fillOpacity = 1,
                   color= "blue",
                   label = paste(expired$trade_name, "|",
                           expired$address_full),
                           group = "Expired Retailers") %>%
  addLayersControl(
    overlayGroups = c("Active Retailers", 
                      "Expired Retailers"),
    options = layersControlOptions(collapsed = TRUE)) %>%
  addLegend(color = "blue", 
            labels = "Expired Retailers",
            group = "Expired Retailers",
            position = "topright") %>%
  addLegend(color = "red", 
            labels = "Active Retailers",
            group = "Active Retailers",
            position = "topright")  %>% 
  hideGroup(c("Expired Retailers")) %>%
  addScaleBar(position = "topleft")

m4 # Print the map

11. Data dictionary

11.1. Data frames generated during updating routines

11.1.1 Pennsylvania

retailers_Socrata_pa - Data loaded via the SODA API for PA

stored_pa - Baseline data set

joined_pa - full join between retailers_Socrata_pa and stored_pa

to_geocode_pa - retailers in need of address geocoding

geocoded_pa - retailers which have gone through the geocoding process

errors_pa - data which might have bad geocodes and/or have resisted geocoding and need to be manually geocoded.

geocoded_pa_fixed - data for which broken geocodes have been manually fixed

11.1.2 New Jersey

retailers_new_nj - Data loaded via excel from NJ gov't

stored_nj - Baseline data set

joined_nj - full join between retailers_new_nj and stored_nj

to_geocode_nj - retailers in need of address geocoding

geocoded_nj - retailers which have gone through the geocoding process

errors_nj - data which might have bad geocodes and/or have resisted geocoding and need to be manually geocoded.

geocoded_nj_fixed - data for which broken geocodes have been manually fixed

11.1.3 Delaware

retailers_Socrata_de - Data loaded via the SODA API for DE

stored_de - Baseline data set

joined_de - full join between retailers_Socrata_de and stored_de

to_geocode_de - retailers in need of address geocoding

geocoded_de - retailers which have gone through the geocoding process

errors_de - data which might have bad geocodes and/or have resisted geocoding and need to be manually geocoded.

geocoded_de_fixed - data for which broken geocodes have been manually fixed

11.2. Canonical names for columns in our data sets

trade_name type: chr example: "1221 INC"

account type: chr example: "2005211534"

expired_y_n type: chr example: "ACTIVE"

expiration_date type: date example: 2020-12-31

lat type: chr example: "39.789164"

lon type: chr example: "-75.697388"

address_full type: chr example: "7441 LANCASTER PIKE HOCKESSIN DE 19707"

county type: chr example: "PHILADELPHIA"

publish_date type: date example: 2018-01-01

license_type type: chr example: "TOBACCO PRODUCTS RETAILER"

state type: chr example: "DE"

12. To do

Next steps

• Update data sets with new NJ data
• Run PA and DE updates
• Incorporate these new data into geoscan routines
• Add NYU routines into our exposure routines - stay events etc
• Begin analyzing existing collections of data
• Add more built environment and risk-terrain functionality